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1. Introduction 

Many inverse scattering problems in imaging are known to be nonlinear. Physically, 
this is a manifestation of the fact that the probing waves do not propagate via well- 
defined trajectories. When such trajectories do exist, the inverse problem can usually 
be linearized as is the case, for example, in single-energy computed x-ray tomography. If 
the probing waves experience scattering and trajectories can not be defined in principle, 
nonlinearity of the inverse problem is practically unavoidable. 

Mathematically, the nonlinearity of the inverse problem can be understood as the 
nonlinear dependence of the measured signal on the quantity of interest. In the case of 
optical tomography (OT), the measured signal is the intensity of light exiting from a 
highly-scattering sample and the the quantities of interest are the absorption and the 
scattering coefficients. The nonlinear nature of the dependence of OT measurements on 
these coefficients is well-known [1,2]. 

Practical approaches to solving nonlinear inverse problems can be divided into two 
broadly defined classes of iterative and analytic methods. Iterative methods, including 
the Newton-type [1,3,4] and Bayesian [5] methods, seek to optimize a certain cost 
function according to an iterative rule which typically requires solving the forward 
problem at each iteration step. The advantage of iterative methods is their generality, 
since they do not require knowledge of the analytical structure of the forward operator. 
Instead, the forward problem is solved at each iteration step numerically. Methods 
of the second class rely on some analytical manipulations with the forward operator. 
This includes various approximate linearization schemes which, generally, work only 
for weak inhomogeneities and methods based on functional series expansions. Thus, 
image reconstruction algorithms based on an inverse scattering series were proposed in 
geophysics (inverse scattering of seismic waves) [6,7], in optical near- field imaging [8], 
and in OT [9]. 

While little is presently known about the convergence of the inverse series, a 
number of results on convergence of the forward series has been obtained. In quantum- 
mechanical scattering theory, Bushell has shown that the Born series converges if the 
potential is too shallow to support at least one bound state [10]. Colton and Kress 
have studied the convergence of the Born series for the scalar wave equation in an 
infinite space [11]. In particular, as part of the proof of Theorem 8.4 of Ref. [11], it 
is shown that the Born series converges if the susceptibility ?7(r) = n^(r) — 1 [n(r) 
being the refractive index] is bounded by \r]{T)\ < 2/{kay, where k = u/c is the wave 
number and a is the radius of the smallest sphere that contains the support of ?7(r). 
We note that Bushell's convergence condition is indirect and, therefore, difficult to use. 
The convergence condition of Colton and Kress is applicable to functions of compact 
support but is not useful at all in the limit /ca — > oo. In addition, it is only applicable 
to scattering by a potential in free space. In this paper, we show that, in the case of the 
diffusion equation used in OT, a simple condition for convergence of the Born series can 
be obtained independently of the medium boundaries. A remarkable property of this 
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condition is that it is also independent of shape or support of the inhomogeneity but 
only depends on its amplitude. Thus, we show that the forward series expansion for the 
Green's function of the diffusion equation in powers of absorptive inhomogeneity (5«(r) 
(the absorption coefficient is decomposed as air) = + 5a(r) where is a constant) 
always converges if 

|5«(r)|<«o. (1) 

A similar condition is obtained for the diffusion coefficient -D(r) = Dq + 5D{v). We 
argue that the independence of this condition on the shape or spatial extent of the 
inhomogeneity is a consequence of the exponential decay of diffuse waves which results 
in weak long-range interactions. This argument will be made more precise in Section [5] 
below and illustrated numerically in Section [71 

The convergence condition ([1]) is obtained independently of restrictions on the 
support of the inhomogeneities or of the nature of medium boundaries. However, if the 
support of the inhomogeneity is contained in a finite ball of radius a and the system is 
embedded in an infinite homogeneous medium, we can repeat the arguments of Ref. [11] 
for the diffusion equation and obtain an even sharper condition on 5a(r). Namely, we 
will show that, for absorbing inhomogeneities and under the conditions stated above, 
the Born series converges if 

1 - (1 + Mexp(-M ^ ' 

where kd = ^Ja^/ Dq is the diffuse wave number (the analog of the wave number k of 
the scalar wave equation). It can be seen that in the limit kda oo, we reproduce 
the condition 6a < Oq, while in the limit k^a — 0, we reproduce Colton and Kress' 
condition 6a < 2ao/{kdaY (note that 6a/aQ is the direct analog of the susceptibility rj 
of the scalar wave equation considered by Colton and Kress). 

The paper is organized as follows. In Section [2] we define the problem of OT, review 
the mathematical formalism that leads to the Born series expansions and introduce the 
relevant notation. In Sections [3] and H] we obtain the convergence condition of the type 
([T]) for absorbing and scattering inhomogeneities, respectively. In Section [5] we generalize 
the Colton and Kress' result for the case of the diffusion equation with an absorbing 
inhomogeneity embedded in an infinite homogeneous medium and derive the convergence 
condition ([2]). In Section [6] we describe a discretization scheme for representation of 
operators by matrices which is used in numerical examples of Section [71 Here the 
analytical results of Section [3l are verified numerically. Finally, Section [HI contains a 
discussion of obtained results. 

Before proceeding with the main content of this paper, we wish to clarify the 
following point. In the text below, we use the terms "multiple scattering of diffuse 
waves" and "interaction" . We are referring to multiple scattering of scalar solutions to 
the diffusion equation from inhomogeneities in its coefficients - not to multiple scattering 
of electromagnetic waves from inhomogeneities in the dielectric susceptibility. The first 
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effect can be viewed as macroscopic, and takes place on much larger scales than the 
second effect. In particular, a macroscopically homogeneous medium with constant 
absorption and diffusion coefficients exhibits no scattering of diffuse waves, although 
the very possibility to describe the electromagnetic energy density by the diffusion 
equation is based on the assumption of strong multiple scattering of electromagnetic 
waves on microscopic physical scales. Similarly, by "interaction" we mean the interaction 
(interference and multiple scattering) of diffuse waves scattered from macroscopic 
inhomogeneities. 

2. Derivation of the Born Series 

Propagation of light in biological tissues is commonly described by the diffusion 
approximation to the radiative transport equation [1,2]. In the case of continuous-wave 
illumination, the following steady-state diffusion equation is used: 

[-V-D(r)V + «(r)]n(r) = g(r) , (3) 

where u is the energy density of the diffuse light inside the medium, q is the source 
function, D = c/[3(/ia + fi'^)], a = cfia and c is the average speed of light in the medium. 
Further, and fi'g are the absorption and reduces scattering coefficients, respectively. 
Reconstruction of the functions /ia(r) and fi'g{r) from a set of boundary measurements 
is the goal of OT. 

Experiments in OT are usually performed with point sources (plane-wave [12] or 
structured [13] illumination have also been proposed). A point source can be written 
as q{r) = qo6{r — r^). Here r^ is the source location on the boundary of the medium. 
A point detector located at can be shown [14] to produce a measurement that is 
proportional to the Green's function of Eq. ([3]), G{rd,rs), which satisfies 

[V ■ D(r)V - a(r)]G(r, r') = -6{r - r') . (4) 

We now decompose a{r) and -D(r) as constant background values Oq, Dq and 
spatially-varying functions Sa{r), SD{r), according to a(r) = + 6a{r) and -D(r) = 
-Do + 6D{r). The background constants are chosen to be equal to the respective values 
of a and D near the medium boundary where these coefficients are either directly 
measurable or known, i.e., by immersing the sample into a matching fluid whose 
optical properties are known. We then obtain the Dyson equation for the Green's 
function [15,16], namely 

G(r, r') = Go(r, r') + J Go(r, r")V{r")G{v" , v')d\ , (5) 

where the integration is over the spatial region occupied by the scattering medium, 
G'o(r, r') is the Green's function for a homogeneous medium with a = aQ and D = Dq, 
i.e., it satisfies 



[DoV'-ao]G'o(r,r') = -(5(r-r') 



(6) 
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and appropriate boundary conditions on the scattering medium boundary, and V{y) is 
given by 



V{v) = V^{v) + Vd{v) , (7) 
V;(r) = -5a{v) , (8) 
VD(r) = -p.5D(r)p. (9) 

Here we have introduced the momentum operator p = — zV. Since p is Hermitian (self- 
adjoint), so is Vd- We note that the Dyson equation ([5]) is vahd for r, r' being inside 
the scattering medium or on its boundary. In the latter case, we can replace r and r' 
by Td and r^. 

In operator notation, the Dyson equation ([5]) is written as 

G = Go + GoVG , (10) 

where V = + Vr, ''^^ the interaction operator. We note that Va is diagonal in the 
position representation and has the matrix elements 

(r|K|r') = -(5a(r)5(r - r') . (11) 



However, Vr, has no position representation. |§| Its matrix elements can be defined in 
the basis of plane waves (in k-space). For example, in an infinite space we can take the 
basis functions to be l^/'k); such that (r|^/;k) = (27r)~^/^ exp(zk ■ r). Then we have the 
following matrix elements (of both Vq, and Vd): 



(^k'|K|V^k) = -5a(k-k') , (12) 
{MyD\i^^.) = -k' ■ k 5D(k - k') , (13) 

where the tilde denotes three-dimensional Fourier transform with respect to the spatial 
variable r. The simple mathematical structure of the above matrix elements suggests 
that the forward and inverse problems are more naturally formulated in k-space, 
especially if the medium boundaries are translationally invariant [14,17]. 

The Born series is obtained by iterating (fTOil starting with G = Gq and has the 

form 

oo 

G = Go + G^VG^ + GqVGoVGq + . . . = Go ^{VGof . (14) 

k=0 

The Born series can also be viewed as the Taylor expansion of the formal solution to 
( |T0|) into a power series in V, 

G={I- GoVy'Go = Go(/ - VGo)-' , (15) 

§ Of course, differential operators in Eq. ([3]) can be approximated by finite differences. However, all 
finite difference sehcmes are non-local (involve several spatial points) and, strictly speaking, can not be 
used to define a position representation of Vd ■ 
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I being the identity operator. 

The derivation of the convergence condition can be obtained directly starting from 
Eq. ( |T5l) . However, a more mathematically elegant approach can be based on an 
analogous formula for the T-matrix. In the T-matrix formalism, one writes the Dyson 
equation (fTOl) as 



G = Go + GoTGo . (16) 
From the identity TGq = VG, we obtain T = VGGq^ or, substituting this into (fT5|) . 

T = V{I -GoVy^ = {I -VGoy^V . (17) 
The Born series for the T-matrix is 



T = V + VGoV + VGoVGqV 



J2^VGo) 



.k=0 



V . 



Note that the series in f|T^ and f|T8|) are identical and, therefore, the convergence 
conditions for the series expansions of G and T are also identical. 



3. The Convergence Condition for Absorbing Inhomogeneities 

The diffusion approximation is valid when /i'^ ^ fia. If, in addition, /i'^ is constant 
inside the sample, then D{r) is also, approximately, constant. This case is practically 
important when the contrast mechanism is directly related to absorption, but not to 
scattering, for instance, in imaging of blood oxygenation levels [18]. 

In this Section, we specialize to the case 6D = 0, Sa 0, so that V = V^- We say 
that the function 8a{v) is physically allowable if 6a{r) > —ao. In the opposite case, the 
total absorption coefficient a(r) = ao + 6a{r) can become negative, which physically 
corresponds to an amplifying medium. 

The derivations presented below are based on the assumption that for any physically 
allowable 6a, the diffusion equation ([3]) has a solution. We also use the fact that if Sa is 
physically allowable and satisfies 6a < ao, then —6a is also physically allowable. While 
we assume on physical grounds that Eq. ([3]) has a solution for every physically allowable 
6a, it can not be stated that if 6a is not physically allowable, then ([3]) has no solutions. 
In fact, ([3]) can have a steady-state solution even if the medium is amplifying in some 
finite spatial region, as long as there also exists a sufficiently strong energy sink 0. For 
this reason, the convergence conditions derived in Sections I31H1 and are sufficient but not 
necessary. 

% li D = Do = const, the diffusion equation ([3]) is mathematically equivalent to the Schroedinger 
equation for a single particle of mass m in the potential U{r) ~ {fi^ /2m)a{r) / Dq- From the analysis 
presented below, it will be clear that the solution to ^ ceases to exists if the potential U{r) is deep 
enough to support at least one bound state. See Ref. [10] for a similar argument. 



On Convergence of Born Series 



7 



3.1. Sign-Definite 6a. 

We start with the simple case of a sign-definite function 5a; (r). Namely, we assume that 
5a{v) does not change sign within its domain (but can be zero). We also assume that 
(5a (r) has no singularities. Then we can write 

V = -aSS , (19) 

where a = ±1 and is a non-negative definite operator, diagonal in the position 
representation. The values of a are a = +1 if (5a > and cr = — 1 if 5a < 0. Then, with 
little algebraic manipulation, we obtain 

T = -aS{I + aSGoSy^S = -(tS(I + aWy^S . (20) 

In the above formula, W = SGqS. The matrix elements of W are given by 

(r|iy |r') = ^AM^\ Go(r, r') ^AM^\ ■ (21) 

The operator W can be viewed as a functional of 6a. We note the following obvious 
property: Vr[75a] = |7|Vr[|5a|], where 7 is a constant. 

W is real and symmetric so that all of its eigenvalues are real. The Born series 
f|T8|) converges if all eigenvalues satisfy < 1 and diverges otherwise. We note that 
the index fi that labels the eigenvalues may not be countable, i.e., if the spectrum of W 
is continuous. Of course, the eigenvalues are not computable analytically in general 
and the above condition is of little practical use. However, we will employ the following 
lemma to obtain conditions on 6a itself: 

Lemma 1 For any physically allowable 6a that does not change sign, (Tu;^[(5a] 7^—1 
for all indices 11. 

Proof For any physically allowable 6a, there is a solution to the diffusion equation ([3]) 
and, correspondingly, a T- matrix. For the T- matrix to exist, the operator / -f- aW in 
(!2T|) must be invertible. But if aw^ = — 1 for at least one eigenvalue, the above operator 
is not invertible. 

In particular, for non- negative functions 6a [a = +1), tf^[5a] 7^—1 and for non- 
positive and physically allowable functions 6a {a = —1), w^[(5a] 7^+1. If the sign 
of a physically allowable 6a is reversed and —6a is still physically allowable, then 
u;^[5a] 7^ ±1. This property holds for all physically allowable functions 6a such that 
6a < aQ. 

We can now state two simple results that set bounds on the spectrum of W. 

Proposition 1 For any physically allowable 6a that does not change sign, Vr[5a] has 
no negative eigenvalues. 
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Proof Let have an eigenvalue w < 0. Choose 7 = Then Ty[7|(5a|] has an 

eigenvalue —1. Since 7|(5a| is non-negative, this is not possible by Lemma 1. 

Proposition 2 If, in addition to the conditions of Proposition 1, 5a < ao, then all 
eigenvalues of iy[(5Q!] are less than unity. 

Proof Let have an eigenvalue w > 1. Choose 7 = —1/w. Then iy[7|(5a|] has an 

eigenvalue +1. Since 7|(5a| is physically allowable and non-positive, this is not possible 
by Lemma 1. 

To summarize, we have found that all eigenvalues of the matrix W = SGqS lie in the 
open interval [0, 1) for all physically allowable functions 5a that satisfy the conditions of 
Proposition 2. Since cr = ±1, we immediately conclude that, under the same conditions, 
the expansion of fl20l) into a power series in W converges. It is further straightforward 
to see that this expansion is identical to ( fTSj) or ( fT4l) . Therefore, we have established 
the following condition for convergence of the Born series: 

Theorem The Born series for the T-matrix or the Green's function converges if (i) 5a 
is physically allowable, (ii) does not change sign inside its domain, and (iii) satisfies 
5a{r) < aQ. 

A remarkable feature of the above condition is that it depends only on the 
upper bound for 5a, but not on its shape or spatial extent. Thus, for example, let 
(5a (r) = A < ao inside some region Q. The Born series will converge independently 
of the shape or linear dimensions of this region. Physically, this can be understood 
by considering the fact that Go{r,r') decays exponentially with the distance between r 
and r'. Therefore, multiple scattering of diffuse waves on large scales is exponentially 
suppressed. Instead, scattering is strong at small scales, when Go{r,r') oc l/|r — r'|. It 
is this short-range interaction that may result in a substantially nonlinear dependence 
of G{V) or T(y) on V. If 5a/ ao is sufficiently large, even locally, the nonlinearity may 
become so strong that the power series expansion of T{V) does not converge. However, 
we have established that this expansion always converges if 5a < aQ. 

We conclude this subsection with the following remark. Proposition 1 is stronger 
than is needed for the derivation of the above convergence condition. The inequality 
> —1 would be sufficient. In fact, we will see below that Proposition 1 holds only for 
operators W whose trace is infinite. If we perform a discretization as is explained 
in Section [6l W becomes a finite-size matrix of zero trace. The scaling property 
iy[7(5a] = |7|H^[|(5a|] does not hold for such matrices. Consequently, some of their 
eigenvalues are negative. However, they are all greater than —1. The proof of this 
statement is very similar to the proof of Proposition 2 and is omitted; instead, we will 
illustrate this fact with numerical examples. 
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3.2. Sign- Indefinite 6a. 

We will now show that the convergence condition formulated in the previous subsection 
holds even if (5a (r) can change sign. 

Before proceeding with the proof, we set the stage for the numerical verification of 
this statement in Section [7] below. Since 5a is now allowed to change sign, we can no 
longer write V = —aSS where a = ±1 and S is real and non-negative definite. Instead, 
we can write, for example, V = —ScSc, where Sc is complex. Analogously to ( l20l) . we 
have 

T = -S,{I + S,GoS,)-'S, = -S,{I + W,)-'S, , (22) 

where Wc = ScGqSc- The matrix elements of Wc are 

{r\Wy) = ^/6^ Go(r, r') y/6a{r') . (23) 

Note that Wc does not depend on the choice of the square root branch in the above 
formula, as long as the same branch is chosen in both square roots. 

Since Wc is complex symmetric and hence non-Hermitian, its eigenvalues are in 
general complex. Therefore, placing bounds on the eigenvalues of Wc is problematic. 
Indeed, the analog of Lemma 1 for Eq. (122|) is 7^ — 1. But this inequality can be 
satisfied trivially if has an imaginary part. Therefore, Eq. (122!) is not useful for 
the derivation of a convergence condition. Instead, we will study eigenvalues of Wc 
numerically in Section [71 Here we will use a different representation for the T-matrix. 
Namely, we can write V = —ST,S where S is still real and non-negative definite but S 
is now an operator rather than a number: 

^^"^"^^=^(^-^)| -1, if Mr)<0. (''^ 

Thus, we can refer to S as the sign operator. Note that E and S commute. After 
straightforward algebraic manipulation, we obtain 

T = -s{^ + SGoSy^s = -s{j: + wy^s . (25) 

In the above equation, is defined by ( 12T1) of Section ISTTj but its domain has been 

generalized to include functions 6a that can change sign. Still, since VTf^o;] = 
and from the results of previous subsection, we know that the eigenvalues of W lie 
in the interval [0,1), as long as \6a\ < a^. Therefore, < 1, where || ■ || is the 

operator norm defined here as = sup[('?/'|W^|'?/')/(^/'|V^)]. On the other hand, from 

the obvious relation = J, we find that = 1. We then write 

(s + wy^ = [s(/ + sw^)]~^ = (/ + ^wy^j: . (26) 

The Born series is obtained by expanding 

oo 

(/ + Siy)-^ = ^(-Eiy)*= . (27) 

fc=0 
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From the operator norm inequality ||y4i?||p < \\A\\p ■ we immediately obtain 

1 1 Ely 1 1 < 1, which is a sufficient condition for convergence of the series fl27j) . This 
completes the proof that the convergence condition of the previous subsection applies 
to functions Sa{r) that can change sign. 



4. The Convergence Condition for Scattering Inhomogeneities 

If fia = const while fi'^ varies, the system is characterized by a scattering inhomogeneity. 
We then have 6a = 0, 6D ^ 0. Obviously, the physically allowable values of 6D satisfy 
SD > —-Do- However, the physical interpretation of what happens if we do allow D{r) 
to become negative is somewhat different. If the source function of Eq. ([3]) is zero in 
the spatial region where D is negative, than the interpretation is that the medium in 
that region is amplifying, similar to the case of absorbing inhomogeneities. But if D is 
negative in a region where the source is nonzero, then, in addition to having amplifying 
medium, the source of energy is turned into a sink. 

We now restrict consideration to a physically allowable SD and state that the 
convergence condition of Section [3] applies to scattering inhomogeneities with the 
substitution — ^ Dq and 6a 6D. The proof of this statement is analogous to 
the proof given in Section [3] and will be only briefly sketched. 

For a general physically allowable 5D, the interaction operator can be written as 
V = Vd = — p ■ SUSp and the symmetric expression for the T-matrix, analogous to 
is 

T = -p- S[J: + SpGop- Sp . (28) 

The operator W = SpGop ■ S is complex but Hermitian, so that all of its eigenvalues 
are strictly real. By considering the special cases of sign-definite SD when E = ±7, 
we obtain bounds on the eigenvalues of W in complete analogy with Section 13. 1[ More 
specifically, the eigenvalues of W all lie in the open interval [0, 1), as long as 6D < Dq. 
We then find that the operator norm of W is less than unity while it is exactly unity 
for E, and, consequently, expansion of fl28l) into a power series converges. 



5. Generalization of Colton and Kress' Result 



Further insight into the convergence properties of the Born series and the strength of 
nonlinearity can be gained by considering the argument similar to the one used by 
Colton and Kress in the proof of Theorem 8.4 of Ref. [11]. The argument is based on a 
direct estimation of the norm ||VGo||oo of the operator VGq that appears in the series 
iHM or (ITSl) . The (necessary and sufficient) convergence condition for the Born series 
is llVGolloo < 1- Of course, estimation of this norm is possible only if Go is known 
analytically. For a medium with boundaries. Go can only be computed numerically. 
Therefore, we will consider below the simple case of free space, so that 
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where = yao/-Do is the diffuse wave number. However, note that the inffuence of 
boundaries can be exponentially small, as is discussed in Section [6] below. 

Next, we specialize to the case of absorbing inhomogeneities, V = Va, where Va is 
defined by ([8]). Assuming that 6a{r) = if r is outside of a sphere of radius a, we have 



iV^Golloo < sup (|5a(r)|) sup (|J(r)|) , (30) 

|r|<a |r|<a 



where 



/(r) = [ GF{r,r')£r' . (31) 

Jr'<a 

The above integral can be easily evaluated to yield 

(32) 



Kr) ^ 



, , ^ . , . exp(/crfr) -exp(-/cdr) 
^1 + kda) exp(-fcda) 



2kdr 

Obviously, the maximum of the above function is at the center of the ball, so that 

sup(|/(r)|) = 7^/(M , /(x) = l-(l+x)exp(-x) . (33) 

|r|<a ^Of^d 

We then immediately arrive at the (sufficient) convergence condition ([2]) for 6a. 

We now examine the two limiting cases k^a — > and kaa — > oo. In the first case, 
we use f{x) ~ for small x and recover Colton and Kress' convergence condition 
6a < 2aQ / {kda^ . In the second case, the domain of 6a is not restricted and we recover 
the result of Section [3|, namely, 6a < aQ with the only difference that we now have a 
strict inequality. The independence of the latter result on /c^a is specific to the diffusion 
equation and results from the exponential decay of diffuse waves. Indeed, we have 
limfc^a^oo fikdo) = 1. However, if we perform the analytic continuation kd — > ik, the 
corresponding limit is limfca^oo \f{'>'ka)\ = ka and the convergence condition becomes 
rj < 1/ka (we have replaced here 6a/ ao by its counterpart 77). This fact illustrates 
the crucial difference in convergence properties of the Born series for propagating and 
diffuse waves. 



6. Discretization 



In any numerical simulations, the operators Go, V must be discretized and 
truncated using some appropriate basis. Here we restrict our attention to absorptive 
inhomogeneities so that V = Va and use the basis of cubic voxels. We note that the 
same discretization method can not be applied to Vd because, as was mentioned in 
Section [2], Vd has no position representation. 

The discretization method described below is analogous to the so-called discrete- 
dipole approximation [19] that has been widely used in electromagnetic scattering by 
nonspherical particles [20,21]. We seek to discretize the integral equation ([5]) in a 
basis of cubic voxels. Instead of working directly with ([5]), it is more convenient 
to first write the Lippmann-Schwinger integral equation for the field u itself. Let 
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u{y) = J G{r,r')q{r')d^r' and M™'^(r) = J Go{r,r')q(r')d^r' . Here M™'^(r) is the incident 
field, i.e., the field that would exist in the absence of inhomogeneities. Using V = V^, 
we obtain the following integral equation for u{r): 

u{r) = u'""(r) - j Go(r, r')5a(r')u(r')ciV . (34) 

We then break up the sample into cubes C„ of side h, volume v = h^, and denote the 
center of each cube by r„. The field ^(r) is approximated by a set of discrete values 
Un = u{rn)- Setting r = r„ in flM|) and representing the volume integral as a sum of 
integrals over each voxel, we obtain 

= -Y, f ^o(r"' r)5a(r)M(r)rfV , (35) 

where u™'^ = M™'^(r„). The above equation is, so far, exact. We now introduce several 
approximations. First, we replace Sa{r)u{r) in the integrand of Eq. ( l35i) by 6amUm, 
where 5a„ = 6a{rn). Second, in all terms with m 7^ n, we replace Go(i"n,r) by 
G'o(r„,rm). We then have 



Un = u'^"" - ^ G'o(r„, Tm)vSamUm. - Q„(5a„Mn , (36) 

Qn= I Go(r„,r)dV. (37) 

Note that the term with m = n has been treated separately because the homogeneous 
medium Green's function Go(r, r') has a singularity at r = r'. The singularity is 
integrable and the quantity Q„ is well defined. However, the computation of Q„ is 
complicated due to the following two factors. First, Gq{y,y') depends on the shape 
of boundaries and on the extrapolation distance £ in a complicated way and is not 
computable analytically in general. Second, the integration in ( 1371) is over a cubic 
volume, while the asymptotic limr^r'[Go(r, r')] oc l/|r — r'| has spherical symmetry. 
The first difficulty is resolved by noting that Gq is a sum of the Green's function in an 
infinite homogeneous space Gp and a contribution due to the boundaries Gb'- 

Go(r,r') = G^(r,r') + GB(r,r') , (38) 

where Gi? is given by (|29|) . Accordingly, we can write Qn as a sum of two contributions, 
Qf and Qbu- Note that the is independent of the index n because the Green's 
function in an infinite homogeneous space is translationally invariant. The term Qbu 
can depend on n because boundaries break translational invariance, so that the integral 
in fl37|) can depend on r„. However, we will argue that Qbu is a small correction to Qp- 
Indeed, GBi^n, r) can be written as a surface integral taken over the medium boundaries 
and has no singularity at r = r„. We estimate that Qbu/Qf ~ {h/ Ln) exp{—kdLn), 
where is the characteristic distance from the point r„ to the medium boundary. We 
assume that all inhomogeneities are localized in a spatial region which is sufficiently far 
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from the medium boundaries. Then the ratio Qbu/Qf is at least of the order of h/Ln, 
if, in addition, k^L ^ 1, this ratio is exponentially small. Therefore, we will neglect the 
term Qbu- The second difficulty is resolved by replacing the integration over the cube 
Cn by integration over a sphere of equivalent volume centered at r„. The radius of this 
sphere is -Req = {S/AnY^^h. With these two approximations, and using (13T1) . (l32l) . we 
have 

1 

where f{x) is defined in ( 1351) . Note that for small x, Qp ^ R^J2Dq. 

Having computed Qp, we can write a self-consistent "coupled-dipole equation" 
which is a discrete approximation to the integral equation (!35ll . § We define "dipole 
moments" dn = —vUnSan, and, after some rearrangement of fl36p . obtain 



Qn = QF = ■^jTrfC'dRccd , (39) 



G'o(r„, Tm)d„ 



(40) 



X» = -T^^ . (41) 

In the above equation, Xn plays the role of polarizability of the ra-th dipole. In the 
absence of interaction, dn = XnU^n^- Note that the polarizability depends on 5an 
nonlinearly due to the presence of the term Qp^otn in the denominator. A nonzero value 
oi Qf can be viewed as a result of interaction of the n-th dipole with itself and therefore 
can be referred to as the dipole self-energy. The physical effect of self-interaction is to 
limit the polarizability. Thus, the maximum (in absolute value) polarizability obtained 
in the limit 5a„ ^ oo is —v/Qp. We note that in the limit kdRcq — ^ 0, QpSan ^ 1. In 
practice, the term QpSan can be small but not zero and should be accounted for. 

We now return to operator notation. Let \d) be an iV-dimensional vector of dipole 
moments dn, n = 1, . . . , N, where N is the total number of voxels. Similarly, we define 
the iV-dimensional vector We then have 

\d) = V4W^') + Gr\d)] . (42) 
Here Va and Gq^ are x A^-matrices with elements 



{n\Va\m) = Xn^nm , (43) 

(n|G7|m) = (1 - 5„„)Go(r„, r^) . (44) 

In the above formula, the superscript "VV" is an abbreviation for "volume-to-volume" 
and is used to emphasize that r„ and are inside the discretized region. The formal 
solution to ( l42i) is 

\d) = {I - V^G^yW^ . (45) 

In the case of the scalar field u(r), a more appropriate term is "couplcd-monopole equation" since the 
quantities dn are, in fact, monopoles. We. however, adhere to the terminology used in electromagnetic 
scattering theory. 
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If there are A'^, discrete sources located at the points r^fc (fc = 1, . . . , N^) and A^^ discrete 

detectors at points r^; (/ = 1, . . . , A''^), we can write within the same precision as was 
used to discretize Eq. ( l35i) : 

= + G^^ (J - KG^) . (46) 
where the matrices G^^, Gq^, Gq^ and G^^ have the following elements: 

(/|G°S|A;) = G(r,,,r,fc) , (47) 

{l\G^^\k) = Go{v,i,v,k) , (48) 

(/|G°V|n) =Go(r,,,r„) , (49) 

(n|Gf |A;) = Go(r„,r,fc) . (50) 



Thus, G°^ and Gg ^ are matrices of size Nd x Ng, Gg^ is of size Nd x N and Gg^ is of 
the size x A^^- The superscripts "VS" and "DV" stand for "source-to-volume" and 
"volume-to-detector" , respectively. 

Eq. P6l) is a discrete approximation to f|T6|) . We can identify 

T = (/-\4Gr)~V« (51) 

as the discrete approximation to the T-matrix while Va and G^^ as discrete A^- 
dimensional approximations to the operators Va and Gq that were considered in 
Sections [2|3l We can further define the square root of Va- For example, if 6an are sign- 
definite, we write Va = —aSS, where 5* is a diagonal matrix with the elements 
Then the T-matrix is written in the symmetric form (l20l) with W = SG^^^ S. In the 
case of sign-indefinite we write the T-matrix in the form (l22l) with Wc = ScG'^^Sc 
and Va = —ScSc (see Section [3721) . 

The T-matrix can be computed by direct inversion of / — V^G^^. This problem 
is well posed and has computational complexity 0(A^^). It should be stressed that 
computation of the T-matrix is completely independent of the sources and detectors and 
only requires knowledge of 6a{r) and the unperturbed Green's function Go(r, r'). Once 
the T-matrix is found, the signal for any source-detector arrangement can be computed 
using (l46i) by direct matrix multiplication, an operation that can be performed with 
computational complexity 0[A^^ min(A'£;, A'^5)-|-A^A'^rfA'^5]. In a situation when the number 
of measurements is approximately equal to the number of unknowns, e.g., A^ ~ NsNd, 
the complexity of matrix multiplication is negligible compared to the complexity of 
computing the T-matrix. 

The T-matrix approach to solving the forward problem has several advantages 
compared to finite differences or finite elements methods. First, only the spatial regions 
where inhomogeneities are supported need to be discretized. In this sense, the method 
is somewhat analogous to methods involving adaptive mesh generation. Second, once 
the T-matrix is computed, the measurable signal can be easily found for an arbitrary 
configuration of sources and detectors. However, unlike the finite difference and finite 
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elements methods, the T-matrix method requires knowledge of Go(r, r') which satisfies 
the proper boundary conditions. We note that Gq can be found analytically for simple 
geometries or, in more complex cases, it can be computed numerically once, e.g., by 
finite differences or the finite-element method. 

We conclude this section by noting that the discretized matrices W and Wc have 
zero trace, unlike their continuous counterparts whose traces are infinite. This is 
due to the renormalization procedure that was employed to remove the singularity of 
Gq{jc, r'). Correspondingly, the sum of all eigenvalues of W or Wc is zero. Some of the 
eigenvalues of W are necessarily negative. In practice, we will see that W has many 
negative eigenvalues of very small absolute value and a much smaller number of positive 
eigenvalues. When < ao, all eigenvalues are located in the unit circle. 

7. Numerical Examples 

We now illustrate the theoretical results of Section [3] with numerical examples using 
the discretization scheme of Section [6l All simulations have been performed in an 
infinite space, so that G'o(r, r') = G^(r, r'), where Gp is given by fl29|) . Physically, 
this corresponds to sources, detectors and the sample being immersed into an infinite 
homogeneous scattering medium. However, even if the sources and detectors are placed 
on the boundary (a diffuse- nondiffuse interface), the replacement of Gq by Gp can 
be a reasonably accurate approximation if the boundaries are sufficiently far from the 
discretized region. Indeed, as was discussed in Section |6l Gq can be written as a sum of 
Gp and Gb-, where the boundary contribution Gb has no singularities when both of its 
argument are inside the medium but not on the medium boundary. Because (^^(r, r') has 
a singularity at r = r', it dominates at small scales. Since the large-scale interaction 
is suppressed due to the exponential decay of diffuse waves, the input of boundaries is 
relatively insignificant for the computation of the T-matrix. However, computation of 
the data function (the measurable signal) according to the formula ( H6|) can depend on 
boundary conditions very strongly. This is because elements of the matrices G^^ and 
G^^ are the Green's functions Goivd, r) and Go(r, r^) where and are located on the 
medium boundary. 

For the specific choice Go = Gp-, the T-matrix depends parametrically on /c^ = 
ao/Do but not on ao and Dq separately. The same is true for W and Wc. The quantity 
kd is known as the diffuse wave number and = 2TT/kd as the diffuse wavelength; 
it gives the inverse scale on which diffuse waves exponentially decay. In all numerical 
examples shown below, sets the physical scale of the problem. The discretization step 
h is not a physical scale; it merely characterizes the precision to which we approximate 
the continuous field u{y) by a set of discrete values Un- 

In the numerical simulations shown below, we have used LAPACK subroutines 
implemented in Intel's MKL library. In particular, we have used the routines DSYEVD 
and ZGEEV for diagonalization of real matrices W and complex symmetric matrices W^ 
respectively. The computation time (on an 4x1.6 GHz Itanium-II HP rx4640 server) 
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Figure 1. Eigenvalues w„, in descending order, vs the eigenvalue number n, for an 
absorbing inhomogeneity of cubic shape of size H — Xd/'2 and various levels of contrast, 
K. The target is discretizcd by 10'^ cubic voxels of size h = Ad/20. 

scaled approximately as 0.5(iV/1000)3sec for SYEVD and 12(A^/1000)3sec for ZGEEV. 
We have also employed the Rayleigh quotient to compute the maximum eigenvalue of 
the real matrix W. This method is quite reliable and can be used to find the maximum 
eigenvalue of matrices with ~ 70, 000 in approximately one minute (once the matrix 
Gq^ is computed, which can take several additional minutes). 

Although we show no directly relevant data, it is interesting to comment on the 
efficiency of computing the T-matrix by direct inversion of the matrix A = I — VaG^^ 
according to (^5T\i . In the case of sign-definite 6a, factorization and subsequent 
inversion of A by the routines DPOTRF and DPOTRI is performed in approximately 
0.14(A^/1000)3sec. For sign-indefinite 6a, the routines DGETRF and DGETRI were 
employed with a computational time of 0.19(A^/1000)^sec. Thus, computation of the 
T-matrix may be a highly efficient method of solving the forward problem of OT and 
can be applicable for discretization involving up to ~ 10^ voxels. We stress that only the 
spatial regions that support inhomogeneities must be discretized. The computational 
disadvantage of the T-matrix approach is that the matrices and A are dense and 
require large storage and fast access to memory. 

7.1. Sign-Definite Case 

We start with the case when 6a{r) does not change sign. Namely, we compute the real 
symmetric matrix W and find its eigenvalues for several shapes of 6a{r). 

The first example is an absorbing inhomogeneity ( "target" ) which has the shape of 
a single cube with side H = Xd/2. It was assumed that 6a{r) = naQ inside the cube and 
is zero outside. The target was approximated by 10^ cubic voxels of volume h^. For this 




Figure 2. Eigenvalues Wm in descending order, vs the relative eigenvalue number, 
n/N, where N is the size of the T-matrix, for an absorbing inhomogeneity of cubic 
shape, contrast k = 1, and various side length H. The discretization step is ft, = Ad/20. 
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Figure 3. Maximum eigenvalue of W, Wmax, for a cubic target of contrast k = 1 
as a function the cube size H (relative to the diffuse wavelength Xd) for different 
discretization. 



discretization, h = Arf/20, kdReq = 0.195 and Qf^o = 0.053. The contrast k was varied 
from 1 to 4. The eigenvalues of W are shown in Fig. [H Note that for the minimum 
physically allowable contrast k = —1, the eigenvalues differ from the case k = 1 only 
very slightly due to the reversal of sign of the term Qp5an in the denominator of (HTl) 
(data not shown). It can be seen that all eigenvalues satisfy ty„ < 1 for k = 1 with a 
large margin. Obviously, the eigenvalues are even smaller for k < 1. 
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Next, we fix the contrast at k = 1 and study the dependence of eigenvalues on the 
size of the cubic target, H. In Fig. [21 we plot eigenvalues for cubes of varying sizes H 
while the discretization step is fixed at /i = Ad/20. It can be seen that the maximum 
eigenvalue Wmax (the one with the lowest relative number) increases with the cube size 
but, for the set of parameters used, does not exceed unity. To study the behavior of 
li'max in a broader range of parameters, we have used the Rayleigh quotient for various 
cube sizes and three different voxel sizes (see Fig. [3]). The Rayleigh method is well suited 
for computing Wmax because of the large gap between the first two eigenvalues. The size 
of the cube was limited (depending on discretization) by the computational restriction 
on N . The maximum value of N used was = 74, 088. Approximately one fourth of all 
data points were verified by full diagonalization, with very good agreement. It follows 
from Fig. [3] that Wmax does not exceed unity for a very broad range of parameters. The 
curves Wma.x{H/ \d) approach unity from below but appear to be unlikely to cross it. Note 
that inhomogeneities of sizes significantly larger than those used in Fig. |3] are rarely, if 
ever, encountered in OT experiments since the typical value of in biological tissues is 
5cm. The visible difference between curves with h = A^/IO and h = Ad/20 is due to the 
presence of the /i-dependent self-energy QpSan in the denominator of ( 14T|) . This term 
is comparable to unity for h = A^/IO but is already small for h = Ad/20. Therefore, the 
difference between the h = Ad/40 and the h = Ad/20 curves is insignificant. Note that 
we expect that discretization with h = Ad/ 10 is too rough to produce accurate results. 
However, the difference (or the lack of it) between the curves Wma.x{H / X^) with different 
h/Xd can not be used per se to verify convergence of the T-matrix with h. 

Since we have performed numerical simulations in infinite space, it is possible to 
compare w^na^iH/ Ad) with the result that can be inferred from the convergence condition 
([2]). To this end, we note the following. The data for Fig. [3] were computed for a cube 
of contrast k = 1. If we increase the contrast by the factor 7, the Born series will 
still converge as long as 7Wmax < 1, or, equivalently, 6a/ ao < 1/wmax- On the other 
hand, the convergence condition ([2]) has the form 6a/aQ < l/f^kda), where /(x) is 
defined by (l33l) and a is the radius of the smallest sphere that circumscribes the cube of 
side H, namely, a = \piHj2. For these two conditions to be consistent, we must have 
Wraiai{H / Xd) < f {iTy/SH / Xd) ■ The latter function is shown as a dotted line in Fig. [3l 

Next, we consider the effects of multiple scattering of diffuse waves between two 
spatially separated absorbing inhomogeneities. To this end, we plot the spectrum of 
eigenvalues of W for two equivalent cubic targets of contrast k = 1 and side H = Ad/2, 
placed side-by-side and separated by the surface-to-surface distance AH. The targets 
were discretized using h = Ad/20, so that each cube was approximated by 10'^ voxels. 
The results are shown in Fig. HI When the cubes are sufficiently far apart {AH = H), 
the interaction is weak and each eigenstate is doubly degenerate (this is in addition to 
the triple degeneracy of some eigenvalues which is due to the cubic symmetry). When 
the cubes approach, the degeneracy is broken by interaction. However, the effect of 
interaction is weak even when the two cubes approach each other very closely. At 
AH = 0, the two cubes merge and form a single parallelepiped. At this point, the 




Figure 4. All eigenvalues Wn of the matrix W (in descending order) vs the eigenvalue 
number n for an absorbing inhomogeneity of contrast k, — 1 in the shape of two 
equivalent cubes of side H = Xd/2 placed side- by-side and separated by the surface- 
to-surface distance AH. Each cube was discretized using h = A(i/20 (10^ voxels per 
cube). 



maximum eigenvalue is increased only by 17% compared to the noninteracting limit. 
The weak interaction of spatially separated inhomogeneities is consistent with the idea 
of exponentially suppressed long-range interaction which was discussed in Section I3.1[ 

7.2. Sign- Indefinite Case 

We now turn to the case of sign-indefinite 6a{r). In this section, we will study the 
complex eigenvalues of the matrix Wc defined in Section I3.2[ We note that, unlike in 
the case of W which is independent of the sign of 6a, Wc[—Sa] = —Wc[Sa]. Note that 
the eigenvalues of Wc change sign when the sign of 6a is inverted. 

The first example considered here is two cubic inhomogeneities similar to those used 
to compute the data points for Fig. HI but now one of them has the negative contrast 
K = — 1. In Fig. El all eigenvalues of Wc for this system are shown as dots in the complex 
plane. When the cubes are sufficiently far apart, the imaginary parts of the eigenvalues 
are very small (~ 10~^ for AH = H). This corresponds to the non-interacting limit, 
when the interaction operator Wc is, approximately, block-diagonal, where each block 
is real symmetric. As the cubes approach, some of the eigenvalues acquire imaginary 
parts. The eigenstates with complex eigenvalues are "hybridized" , i.e., they are collective 
eigenstates of the two interacting objects rather than "pure" eigenstates of each object 
taken separately. However, the hybridization is weak. Imaginary parts of the eigenvalues 
do not exceed 0.0015 in absolute value. Again, this is in agreement with the idea of 
exponentially-suppressed long-range interactions. 
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Figure 5. All eigenvalues of the matrix Wc for two cubic inhomogeneitics of equal 
sides H = O.SAd, placed side-by-side and separated by the surface-to-surface distance 
AH. One cube has contrast k ~ +1 and the other k = —1. Discretization: h = Ad/20. 



Next, we consider a layered structure of fifteen tliin square layers of thickness h 
and alternating contrast k = ±1 sandwiched on top of each other to form a cube of 
side H = 0.75Arf. The discretization step is still h = A/20. The eigenvalues of Wc 
are shown in Fig. [61 The displayed data indicate that there are hybridized eigenstates 
(those with complex eigenvalues) and eigenstates associated with an isolated thin layer 
and almost unaffected by the interaction (with almost purely real eigenstates). Overall, 
the absolute values of all eigenstates do not exceed 0.05. In this case, the matrix W is 
negligibly small compared to / and can be neglected. This corresponds to the first Born 
approximation, i.e., T = V. Thus, multiple scattering of diffuse waves for this layered 
structure is quite weak and can be neglected with little loss of precision. 

The final example is one cubic inhomogeneity embedded inside another. Namely, 
a cube of size llh x llh x llh with contrast k, = —1 was "coated" by a larger cube 
of size 21h x 21h x 21h with contrast k, = +1. The contrasts in the inner and outer 
cubes were not additive, so that n = —1 in the interior and k = +1 in the exterior 
of the structure. The discretization step was h = Xd/20, so that the outer cube side 
was i^out = l-05A(i; the inner cube side was ifjn = O.SSA^. The eigenvalues of the 
matrix Wc for this structure are shown in Fig. [71 Note that the vertical scale in this 
figure is the same as in Fig. [6l but the horizontal scale is ten times larger. Thus, while 
multiple scattering of diffuse waves inside each component (e.g., within the regions of 
positive or negative contrast) is much stronger than in the case of the layered structure 
of Fig. [6l hybridization is much weaker. The hybridized eigenvalues can be seen near 
the origin of the complex plane and are all very small in magnitude. At the same time. 
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Figure 6. All eigenvalues of the matrix Wc for the layered absorptive inhomogeneity 
described in the text. 
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Figure 7. All eigenvalues of the matrix Wc for the absorptive inhomogeneity in the 
shape of two embedded cubes described in the text. 



the eigenvalues that are relatively large in magnitude are almost purely real, which is 
characteristic for weak interaction between regions with positive and negative contrasts. 

8. Discussion 

In this paper, we have derived a sufficient condition for convergence of the Born series 
for the forward operator of optical tomography. The condition is quite simple and 
states that the series converge if the relative deviation of the absorption coefficient from 
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its background value 5a(r)/ao does not exceed unity, independently of the support 
of (5a (r). A similar condition was obtained for scattering inhomogeneities which are 
manifested by a spatially inhomogeneous diffusion coefficient. We have considered 
absorbing and scattering inhomogeneities separately; the situation when the absorption 
and the diffusion coefficients can vary in space simultaneously is not discussed in this 
paper. We argue that the convergence condition depends only on the amplitude but not 
on shape of 5a (or 5D) due to the exponential spatial decay of diffuse waves. Because 
of this decay, multiple scattering is suppressed on large scales. We emphasize again that 
we discuss here multiple scattering of diffuse waves - scalar solutions to the diffusion 
equation ([3]) - not electromagnetic multiple scattering which happens at much smaller 
physical scales. In the case when (5a(r) has a compact support in a ball of radius a, a 
sharper convergence condition has been obtained (formula ([2])), which is a generalization 
of the result previously obtained for the scalar wave equation [11]. A crucial difference 
between the convergence condition for propagating and diffuse waves is revealed in the 
limit a — > cxD, as is discussed in Section |5l 

An interesting consequence of the convergence condition is that the nonlinearity of 
the inverse problem of optical tomography can be controlled if the constant ao can be 
controlled. Thus, increasing results in effective linearization of the inverse problem. 
Theoretically, can be chosen arbitrarily. However, the ill-posedness of the linear 
inverse problem tends to increase with a^. This reveals an interplay between the ill- 
posedness of the linearized inverse problem and the degree of nonlinearity of the full 
inverse problem (before linearization). Note that in experiments, ao can be tuned, for 
example, by changing the composition of an index-matching fluid. 

We have performed numerical simulations for absorbing inhomogeneities. All 
numerical data are in agreement with the analytical results of this paper. We have 
found that the derived convergence condition is satisfied for a very broad range of 
parameters which are accessible in numerical experiments. We have also found that the 
effects of multiple scattering between spatially separated inhomogeneities such as two 
separate cubes is quite weak. This is again a consequence of the exponential decay of 
diffuse waves. Interaction of inhomogeneities whose contrasts have different signs was 
found to be especially weak. Thus, for the layered structure discussed in Section \7 .2\ 
the interaction is insignificant and the first Born approximation can be used with high 
accuracy - even though the object is a layered cube of size H = O.TSA^. 

While we have found no substantial interaction between spatially separated 
inhomogeneities, nonlinearity can become strong in bulk inhomogeneities of large 
spatial extent or high contrast. In this case, the nonlinearity results from short-range 
interactions. Here two voxels can strongly interact with each other even if they are far 
apart, provided that there is a continuous path of other voxels connecting them. 

Another aspect of the paper that deserves comment is the independence of the 
results on source-detector orientation. Indeed, it may seem natural that two absorbing 
cubes that block the line of sight will have more effect on the measured signal than 
the same two cubes rotated so that only one of them blocks the line of sight. In fact, 
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convergence or divergence of the Born series can be influenced, to a certain extent, by the 
source-detector arrangement. Indeed, calculation of the measurable signal according to 
( 146|) involves multiplication of the T-matrix by G^^ and G^^ from left and right. These 
matrices are source- and detector-dependent. It can happen that the matrix W has an 
eigenvalue larger than unity so that the Born series for the T-matrix diverges, but the 
corresponding eigenvector has a zero projection on either G^^ or G^^. Then the Born 
expansion of the Green's function G^^ will converge for the selected source-detector 
configuration. However, if the Born series converges for the T-matrix, it converges for 
all possible source-detector pairs. 

Finally, our results pertain only to convergence of the forward series. Analogous 
results on the convergence of the inverse series are not yet known. 
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